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Building  Adaptive  Multiresolution  Schemes 
within  Harten’s  Framework 


Francesc  Arandiga  and  Rosa  Donat 


Abstract.  We  consider  the  cell-average  framework  described  by 
A.  Harten  in  [5],  and  build  the  prediction  operator  using  two  nonlinear 
interpolation  techniques.  We  test  the  resulting  nonlinear,  adaptive,  mul- 
tiresolution scheme,  and  compare  it  with  a linear  scheme  of  the  same 
accuracy.  The  nonlinear  prediction  processes  we  develop  can  also  be  used 
in  the  context  of  iterative  refinement.  Numerical  tests  show  that  this  is 
also  a viable  alternative  for  piecewise  smooth  data. 


§1.  Introduction 

The  goal  of  a multi-scale  decomposition  of  a discrete  set  of  data  is  a “re- 
arrangement” of  its  information  content  in  such  a way  that  the  new  discrete 
representation,  exactly  equivalent  to  the  old  one,  is  more  “manageable”  in 
some  respects.  Some  of  the  best  known  applications  of  multi-scale  decomposi- 
tions derive  from  their  compression  capabilities:  a multiresolution  representa- 
tion of  a function,  i.e.,  of  a discrete  set  which  represents  the  function  in  some 
sense,  can  be  highly  compressed  with  minimal  loss  of  information  content. 
Precisely  because  of  this  potential,  multi-scale  techniques  have  an  emergent 
role  in  numerical  analysis,  where  the  multi-scale  idea  has  been  used  success- 
fully over  the  years,  from  multigrid  techniques  to  hierarchical  bases  in  finite 
element  spaces  or  subdivision  schemes  in  Computer-Aided  Design  (CAD). 

In  the  late  80’s  and  early  90’s,  ideas  from  all  these  fields,  together  with 
a wide  experience  in  the  numerical  solution  of  Hyperbolic  Conservation  Laws 
(HCL)  lead  Ami  Harten  to  develop  a General  Framework  for  multiresolu- 
tion representation  of  discrete  data.  The  building  blocks  of  a multiresolution 
scheme  a la  Harten  are  two  operators  which  connect  discrete  and  continuous 
data:  The  discretization  operator  obtains  discrete  information  from  a given 
signal  (belonging  to  a particular  function  space)  at  a given  resolution  level; 
the  reconstruction  operator  produces  an  approximation  to  that  signal  (in  the 
same  function  space)  from  its  discretized  values. 


Curve  and  Surface  Fitting:  Saint-Malo  1999  19 

Albert  Cohen,  Christophe  Rabut,  and  Larry  L.  Schumaker  (eds.),  pp.  19-26. 

Copyright  ©2000  by  Vanderbilt  University  Press,  Nashville,  TN. 

ISBN  0-8265-1357-3. 

All  rights  of  reproduction  in  any  form  reserved. 


20 


F.  Arandiga  and  R.  Donat 


Harten’s  point  of  view  is  that  the  way  in  which  the  discrete  data  is  gener- 
ated, i.e. , the  discretization  process,  determines  its  nature  and  should  provide 
an  adequate  setting  for  a multiresolution  analysis.  Once  the  setting  is  spec- 
ified, the  choice  of  an  appropriate  reconstruction  operator  provides  the  key 
step  to  the  construction  of  a multiresolution  scheme. 

The  reconstruction  process  lies  at  the  very  heart  of  a multiresolution 
scheme  built  a la  Harten , and  adaptivity  can  be  introduced  in  the  multires- 
olution scheme  at  this  level.  A nonlinear,  adaptive  reconstruction  technique 
which  fits  the  approximation  to  the  local  nature  of  the  data  will  lead  to  a 
nonlinear  adaptive  multiresolution  algorithm  with  improved  compression  ca- 
pabilities. 

The  aim  of  this  study  is  to  examine  a particular  class  of  nonlinear  adaptive 
multiresolution  schemes,  those  using  the  Essentially  Non  Oscillatory  (ENO) 
interpolatory  techniques  of  [6]  in  the  reconstruction  step.  Numerical  experi- 
ments [1,2]  show  that  ENO-MR  schemes  have  larger  compression  rates  than 
linear  ones  when  the  original  signal  or  image  is  composed  of  smooth  parts 
joined  together  by  singularities.  ENO  techniques  can  be  used  to  construct 
very  accurate  interpolants,  which  in  turn  lead  to  multiresolution  schemes  with 
high  compression  capabilities.  When  the  original  signal  is  geometric,  nonlin- 
ear schemes  can  be  used  as  loss-less  compression  techniques,  and  we  show 
some  application  of  this  in  the  last  section  of  this  paper. 

The  nonlinear  prediction  process  can  be  used  also  in  the  context  of  sub- 
division refinement.  This  amounts  to  setting  to  zero  all  scale  coefficients  and 
using  the  prediction  operator  to  proceed  by  dyadic  refinement.  Preliminary 
tests  show  that  these  nonlinear  subdivision  schemes  lead  to  non-oscillatory 
limiting  functions  when  applied  to  piecewise  smooth  data  with  jumps,  and 
open  up  an  interesting  alternative  for  iterative  refinement  of  piecewise  smooth 
data. 


§2.  Cell  Average  Multiresolution  Analysis 

When  dealing  with  discrete  data  coming  from  a piecewise  smooth  function, 
the  simplest  discretization  process,  that  of  considering  the  point-values  of  the 
function,  might  not  be  well  defined,  especially  at  jump  discontinuities.  On 
the  other  hand,  the  discretization  by  cell-averages  procedure  acts  naturally 
on  the  space  of  integrable  functions,  and  it  provides  a more  adequate  setting 
to  deal  with  piecewise  smooth  signals.  Because  of  this,  we  shall  carry  out  our 
numerical  study  within  the  cell-average  framework. 

Images  are  considered  here  as  two-dimensional  signals,  and  we  use  the 
usual  tensor-product  approach  to  design  our  two-dimensional  algorithms. 
Thus,  we  only  describe  the  essential  features  of  the  one-dimensional  setting 
for  the  sake  of  completeness.  The  interested  reader  can  find  the  missing  details 
in  this  section  in  [2]  or  [5]. 

Let  us  consider  a set  of  nested  dyadic  grids  in  [0,1]: 

Xk  = Kfc}^0,  xk=ihk,  hk  = 2-k/No,  Nk  = 2kN0,  k = L,...,  0, 
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where  No  is  some  integer.  The  discretization  by  cell  average  operator  is  defined 
as  follows: 


Vk:Ll[ 0,1]  — > Vfc 


ft  = ( Vkf)i  = ^ f ' f(x)dx,  1 <i<Nk, 


(i) 


where  L^O,  1]  is  the  space  of  absolutely  integrable  functions  in  [0, 1]  and  Vk 
is  the  space  of  sequences  with  Nk  components. 

Due  to  the  relation 


“ = srr  [L ' '<*>* = k if  n*)dx  - + '*>• 


it  is  easy  to  see  that  fc  = L — 1, . . . , 0,  can  be  evaluated  directly  from 

{fi}7=\  without  using  explicitly  (1)  (i.e.,  without  knowledge  of  the  original 
function  /(x)). 

To  define  an  appropriate  reconstruction  operator  for  this  setting  (in  fact, 
a whole  family  of  them),  we  consider  the  sequence  {Fk}  on  the  fc-th  grid 
defined  from  the  cell  values  {ft}  as  follows: 

Fi  = hk  £ ft  = r fWdx  = ^ fi=  (2) 

4=1  J0  hk 

The  function  F(x)  (e  C[0, 1])  is,  in  fact,  a primitive  of  the  original  function 
/(x),  and  the  sequence  {Fk}  represents  a point  value  discretization  of  F(x) 
on  the  fc-th  grid  (with  Fq  — 0).  Notice  that  (2)  establishes  a one-to-one 
correspondence  between  {/*}^0  and 

Let  us  denote  by  T(x;Fk~1)  an  interpolatory  reconstruction  of  the  set 
{Ffe~1}  on  the  grid  Xfc_1,  i.e.,  T(x*-1;  Fk~1)  = Fk~x.  Then,  we  obtain  an 
approximation  , /*,  to  fk  using  (2)  as  follows: 

ft  = (I(x?,Ffe-J)  -!(**_!,  F*-1))/^.  (3) 

Since  Fk{  = F{xki ) = F(xf_1)  = we  obtain 

fti—i  = and  fli  = ±(Fk-1-l(xki_1,Fk-1)). 

(4) 

Let  us  define  the  prediction  errors  as  ek  :=  fk  — fk.  Using  (2)  and  (4),  we 
easily  obtain 

/£-!  = -(/£-/£)• 

Thus,  we  can  simply  store  only  the  prediction  errors  with  odd  indexes;  these 
are  the  scale  coefficients,  dk  = e.ki_l,  of  the  multiresolution  transform. 
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The  multiscale  decomposition  of  the  original  data  fL  is  described  by  the 
encoding  algorithm: 

{Do  k = L, . . . , 1 

JLl  = WLi  + fL  1<  «<**-!  (5) 

dki  = f^i—i - (x(4-r if*'1) - Fj^yht  1 <i<  AT*-!. 

We  recover  the  original  data  with  the  decoding  algorithm: 

{Do  k = 1, . . . , L 

fL  1 = (1(4 _j;  f*-1)  - FLi)/hk  + dk  1 < * < Nk. I (6) 
fli  = 2 ft1  - fLi  1 < i < Nk-1. 

In  our  study  we  consider  only  local  interpolation  techniques  with  La- 
grangian  polynomials,  i.e., 

I(x;Fk)  = qi(x]Fk)  xe[xlvxk), 

where  Qi(x\  Fk)  is  a polynomial  of  degree  r satisfying  qi(xk_1~,  Fk)  = Fk_1  and 
qi(xi\  Fk)  = Fk. 

When  the  stencil  of  points  used  to  construct  qi{x)  is  symmetric  around 
the  ith  interval  (i.e.,  r = 2s- 1,  S = {4S, ' ’ ' ixi+s- 1})>  we  obtain  a centered 
interpolation  technique.  Centered  interpolation  techniques  are  very  often 
used  in  approximation  theory  because  they  minimize  the  interpolation  error, 
thus  leading  to  very  accurate  reconstructions  of  smooth  signals.  It  turns  out 
that  the  multiresolution  schemes  obtained  from  (5)  and  (6)  with  Lagrangian 
piecewise  polynomial  centered  interpolation  techniques  are  equivalent  to  the 
Biorthogonal  Wavelet  (BOW)  schemes  of  [4]  (with  the  box  function  as  the 
scaling  function). 

The  compression  properties  of  BOW  schemes  have  been  widely  analyzed 
in  the  literature,  but  from  an  approximation  theory  standpoint,  it  is  very 
easy  to  study  the  behavior  of  the  coefficients  in  terms  of  the  smoothness  of 
the  underlying  signal  and  the  accuracy  of  the  interpolation  technique.  Notice 
that  the  scale  coefficients  dk  are  related  to  interpolation  errors  at  the  odd 
points  of  the  A:-th  grid.  In  fact, 

4 = (4-i-I(4-i; 

Thus,  if  f(x)  is  sufficiently  smooth  at  [xkZ\,  we  have  dk  = ). 

However,  the  presence  of  an  isolated  singularity  x,t  6 induces  a 

loss  of  accuracy  in  the  polynomial  pieces  whose  stencils  cross  the  singularity. 
The  accuracy  loss  is  related  to  the  strength  of  the  singularity  as  follows  [2]:  if 
[f^Ld  = /^(xd+)  ~ f^pHxd+)  — 0(1)  (p  < r),  and  / is  smooth  everywhere 
else,  we  have 


d 


k 

l ~ 


o{[f{p)])K- 1,  f = « -i, 

0(hrk_1),  otherwise. 


(7) 
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Thus,  centered  interpolation  techniques  lead  to  relatively  large  regions  of 
poor  accuracy  around  singularities,  and  therefore  to  large  detail  coefficients 
at  those  locations  where  the  accuracy  loss  takes  place.  The  consequence  is  a 
loss  in  efficiency  for  the  multiresolution-based  compression  scheme. 

It  seems  reasonable  to  improve  the  efficiency  of  the  multiresolution-based 
compression  scheme  by  improving  the  accuracy  of  the  interpolatory  technique 
used  in  the  reconstruction  step.  Notice  that  when  the  convex  hull  of  the  sten- 
cil used  to  construct  a polynomial  interpolant  is  contained  within  a region  of 
smoothness  of  the  underlying  signal,  the  interpolation  error  (and  the  corre- 
sponding detail  coefficient)  becomes  small.  Thus,  it  is  clear  that  the  key  point 
is  the  construction  of  polynomial  pieces  that  avoid  the  singularity. 

In  the  literature  related  to  the  numerical  solution  of  conservation  laws, 
where  discontinuities  can  spontaneously  develop,  we  find  an  interpolation  pro- 
cedure with  all  the  features  we  need:  Essentially  Non  Oscillatory  (ENO)  in- 
terpolatory techniques  [6]  (it  is  not  surprising  that  Harten  was  one  of  the 
developers  of  these  techniques). 

ENO  interpolatory  techniques  lead  to  piecewise  polynomial  interpolants 
that  are  fully  accurate  except  in  those  intervals  that  contain  singularities. 
The  essential  feature  of  ENO  interpolatory  techniques  is  a stencil  selection 
procedure  that  attempts  to  choose  each  stencil  Si  within  the  same  region  of 
smoothness  of  F(x).  The  stencil  selection  process  uses  the  divided  differences 
of  the  discrete  set  to  be  interpolated  as  smoothness  indicators:  Large  divided 
differences  indicate  a possible  loss  of  smoothness.  The  selection  process  is 
such  that  it  tends  to  look  away  from  large  gradients,  when  this  is  feasible. 

ENO  interpolatory  techniques  are  nonlinear,  because  the  stencil  used 
to  construct  each  polynomial  piece  depends  on  the  function  being  interpo- 
lated. When  the  singularities  are  sufficiently  well  separated  (this  means  that 
there  are  at  least  r + 1 points  in  each  smoothness  region),  ENO  techniques 
lead  to  stencils  such  that  (assuming  the  singularity  is  located  at  the  ith  cell) 
<Sj_i  n <S;+i  = 0.  Hence,  the  detail  coefficients  satisfy 

^ f ' = *.  (8, 

1 | 0(^_i),  otherwise. 

Thus  ENO  interpolants  have  a nearly  optimal  high  accuracy  region,  which 
should  in  turn  improve  the  efficiency  of  the  corresponding  multiresolution- 
based  compression  algorithms. 

The  case  of  a corner  of  / (i.e.,  a jump  in  /')  is  especially  interesting 
because  it  is  possible  to  construct  an  even  better  (in  terms  of  local  accuracy) 
interpolant:  the  ENO-SR  interpolant. 

The  Subcell  Resolution  (SR)  technique  (also  due  originally  to  Harten  [3]) 
allows  us  to  obtain  an  approximation  to  the  location  of  an  isolated  corner  in  a 
continuous  function  up  to  the  order  of  the  truncation  error.  The  approximated 
value  is  then  used  to  modify  locally  (in  the  interval  where  the  discontinuity 
lies)  the  definition  of  the  piecewise  polynomial  interpolant  in  such  a way  that 
the  interpolation  error  is  small  except  for  an  0(hr+1)  band  around  the  corner. 
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Fig.  1.  Left  original,  right  coarse  version. 


Fig.  2.  Left  linear,  right  ENO. 

Recall  that  in  the  cell-average  setting,  the  interpolation  process  is  applied 
to  the  primitive  function.  Since  jumps  in  f(x)  become  corners  in  F(x),  using 
an  ENO-SR  interpolant  in  the  reconstruction  step  lead  to  detail  coefficients 
satisfying 

df  = 0(hrk_1 ),  except  when  x^i-i  ~ xj.  (9) 

The  SR  technique  is,  thus,  appropriate  to  increase  the  efficiency  of  mul- 
tiresolution-based compression  algorithm  for  piecewise  continuous  functions 
with  jumps. 


§3.  Numerical  Experiments 

Let  us  consider  a purely  geometric  image  as  shown  in  Figure  1 (left),  and 
apply  to  it  the  tensor-product  version  of  algorithm  (5).  We  consider  piecewise 
polynomial  interpolants  of  degree  4,  thus  the  accuracy  of  the  reconstruction  in 
the  cell-average  framework  is  3.  In  Figure  2 we  display  the  location  of  non-zero 
scale  coefficients  in  the  multiresolution  representation  of  the  signal.  When 
using  the  ENO-SR  technique,  and  because  all  discontinuities  are  “aligned 
with  the  (tensor-product)  grid” , all  scale  coefficients  are  zero.  This  is  a direct 
consequence  of  the  fact  that  the  ENO-SR  reconstruction  commits  no  error  at 
the  odd  points  in  each  one  of  the  resolution  levels  considered  (L  = 4 in  this 
example).  In  the  case  of  the  ENO  scheme,  the  scale  coefficients  at  the  highest 
resolution  level  are  all  zero.  This  is  a consequence  of  the  nature  of  the  data 
(the  point- values  of  the  signal  at  the  highest  resolution  level),  which  locates 
all  discontinuities  at  the  cell  end-points.  The  ENO  technique  is  perfectly 
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Fig.  4.  Reconstruction  from  the  coarse:  linear,  ENO,  ENO-SR. 


accurate  when  the  discontinuity  is  located  at  a grid  point.  The  technique  is 
fully  accurate  at  all  the  lower  resolution  levels  except  at  the  interval  where 
the  discontinuity  is  located.  Thus,  there  is  only  one  scale  coefficient,  located 
at  the  point  which  is  closest  to  the  singularity.  This  should  be  compared 
with  the  3 scale  coefficients  per  singularity  obtained  with  the  linear  scheme 
(the  lines  showing  the  location  of  non-zero  scale  coefficients  are  thicker  for 
the  linear  scheme).  For  the  sake  of  comparison,  the  number  of  non-zero  scale 
coefficients  in  each  case  is:  Linear  8554,  ENO  1688,  ENO-SR  0. 

We  turn  next  to  the  nonlinear  subdivision  scheme  obtained  by  considering 
the  ENO  and  ENO-SR  in  the  prediction  process.  In  Figure  3 we  show  a 
univariate  process.  Starting  from  the  the  cell-averages  of  a piecewise  smooth 
function  at  a very  coarse  level  (16  points),  we  proceed  by  dyadic  refinement 
until  we  obtain  1024  data.  The  numerical  results  clearly  indicate  that  no 
overshoots  or  undershoots  are  obtained  with  the  non-linear  techniques.  Again, 
the  excellent  properties  of  the  ENO-SR  technique  in  terms  of  approximation 
lead  to  the  best  results. 

In  Figure  4 we  show  a simple  multivariate  test,  the  reconstruction  of  the 
geometric  figure  considered  before  from  a very  coarse  representation  (right  in 
Figure  1).  The  Gibbs-like  oscillations  typical  of  linear  schemes  in  the  presence 
of  discontinuities  lead  to  the  blurring  of  the  edges  observed  in  Figure  4 (left). 
There  is  no  blurring  in  the  reconstructed  image  obtained  with  the  nonlinear 
techniques.  Again,  the  ENO-SR  scheme  leads,  in  this  simple  case,  to  the 
exact  original  image.  One  dimensional  cuts  of  the  reconstructed  figures  are 
displayed  in  Figure  5. 
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Fig.  5.  Horizontal  cuts  linear,  ENO,  ENO-SR. 
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